package scu.maqiang.matlabcases;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class MinimalSurface {

	public static void main(String[] args) {
		String meshFile = "matlabMeshes\\Circular.msh";
		Mesh2T3 mesh2D = new Mesh2T3().readMesh(meshFile, MeshFormat.msh);
		FES2T31 fs = new FES2T31(mesh2D);
		SRMatrix A = new SRMatrix(fs.GetNdof());	
		double[] RHS = new double[fs.GetNdof()];
		double[][] grad = new double[2][fs.GetNdof()];
		double[] coef = new double[fs.GetNdof()];
		Arrays.fill(coef, 1.0);
		double[][] NLCoef = {coef};
		NewIterSSolver solver = new NewIterSSolver();
		double[] x0 = new double[fs.GetNdof()];
		ScalarFunc fixBd = (xy, label, param) -> xy[0] * xy[0];
		ScalarFunc[] bd = {fixBd};
		double[] x = new double[fs.GetNdof()];
		fs.applyBC_MBN(RHS, Direct.All, bd, null, 1, 2, 3, 4);
		while(true) {
			fs.assembleStiff(NLCoef, BVPType.COMMON, A);
			fs.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
			solver.setMatrix(A);
			solver.PCGSSOR(RHS, x, 1.5, 1);
			if(MVO.distance(x0, x) < 1.0e-10) {
				break;
			}
			fs.computeGradient(x, grad);
			Arrays.setAll(x0, i -> x[i]);
			Arrays.setAll(coef, i -> 1.0/Math.sqrt(1 + grad[0][i] * grad[0][i] + grad[1][i] * grad[1][i]));
			A.clear();
		}
		mesh2D.toTecplot("matlabResults\\MinimalSurface.dat", x);
	}

}
